% section 6.4.1 'robustness' / Income Profiles 
% this file shows the differences in CE for GHHW and GLTHI 
% using PSID income and allowing for savings

clear;
load('./results/allstates_main_psid.mat')
load('./results/model_parameters_main.mat')
addpath functions

    MU = model_parameters.MU;
    CAPH = model_parameters.CAPH;
    gamma = model_parameters.gamma;
    rho = model_parameters.rho;
    r = model_parameters.r;
    T = size(MU,3);
    age = [25:(25+T-1)]';

% load objects needed

    Pspot_N = simulation_output.Pspot_N;
    Ppaid_N = simulation_output.Ppaid_N;
    Wio_abi= simulation_output.Wio_abi;
    Wio_rea= simulation_output.Wio_rea;
    Cg_abi = simulation_output.Cg_abi;
    Cg_rea = simulation_output.Cg_rea;
    iX = simulation_output.iX;
    Nstates = size(Pspot_N,1)-1;
    Ppaid_spot = simulation_output.Ppaid_spot;
    nInd = size(iX,2);
    T = 70;
    
% compute results    
    % Ed 13
    % maximum allowed level of assets (for grid) under each contract
    amax_phi = 180;
    amax_st = 800;
    [C_PS_abi,A_PS_abi,CE_PS_abi,C_PS_PHI_abi,A_PS_PHI_abi,CE_PS_PHI_abi] = analysis_PS(Wio_abi,CAPH,MU,Ppaid_N,Pspot_N,iX,Nstates,rho,r,gamma,T,amax_phi,amax_st);
    
    % Ed 10
    amax_phi = 130;
    amax_st = 800;

    [C_PS_rea,A_PS_rea,CE_PS_rea,C_PS_PHI_rea,A_PS_PHI_rea,CE_PS_PHI_rea] = analysis_PS(Wio_rea,CAPH,MU,Ppaid_N,Pspot_N,iX,Nstates,rho,r,gamma,T,amax_phi,amax_st);
    
 % graph (doesn't go to paper)

      Cons_PHI_abi =Wio_abi-Ppaid_N;
      Cons_PHI_abi_mean = sum(Cons_PHI_abi.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2); 
      Cons_PHI_PS_abi_mean = sum(C_PS_PHI_abi.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2)/1000; 
      Cons_HHW_abi =Consumption_HHW(Cg_abi,iX);
      Cons_HHW_abi_mean = sum(Cons_HHW_abi.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2); 
      
      
      plot(age,Cons_PHI_PS_abi_mean,age,Cons_PHI_abi_mean,age,Cons_HHW_abi_mean)
      h = legend('GLTHI savings','GLTHI','GHHW');
      set(h,'Location','northwest');
      set(h,'fontsize',14);
      
      saveas(gcf,'./results/PHI_PHI_savs_HHW_psid.png');
    
% table
    % main results to compare  
    CE_GHHW_abi = simulation_output.CE_HHW_abi;
    CE_FB_abi = simulation_output.maintable.CE_FB(2);

    CE_GHHW_rea = simulation_output.CE_HHW_rea;
    CE_FB_rea = simulation_output.maintable.CE_FB(1);

  
    table_precsavings = array2table([CE_FB_rea CE_PS_rea CE_GHHW_rea CE_PS_PHI_rea ;...
                     CE_FB_abi CE_PS_abi CE_GHHW_abi CE_PS_PHI_abi ]);

    table_precsavings = renamevars(table_precsavings,["Var1","Var2","Var3","Var4"],...
                        ["CE_FB","CE_ST","CE_GHHW","CE_GLTHI"]);

    table_precsavings.gainsGHHW = (table_precsavings.CE_GHHW-table_precsavings.CE_ST)./(table_precsavings.CE_FB-table_precsavings.CE_ST);
    table_precsavings.gainsGLTHI = (table_precsavings.CE_GLTHI-table_precsavings.CE_ST)./(table_precsavings.CE_FB-table_precsavings.CE_ST);
    table_precsavings.metric1 = table_precsavings.gainsGLTHI./table_precsavings.gainsGHHW;
    table_precsavings.metric2 = (table_precsavings.CE_GHHW-table_precsavings.CE_GLTHI)./(table_precsavings.CE_GHHW);

    writetable(table_precsavings,'./results/table_precsavings_psid.xlsx');

                 
                 